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FCREWORD 

In  lats  fall  of  1952,  Dr.  George  B.  Dantzig  at  The  RAND  Corporation 
requested  that  a  piK)grainmer  be  assigned  to  work  with  him  in  developing 
machine  computing  procedures  for  the  simplex  method.  The  writer  was  given 
the  assignment,  which  has  since  grown  into  a  full-time  job,  and  has  been 
assisted,  since  late  1953,  by  Miss  Leola  Cutler  and  others  of  RAND's  pro¬ 
gramming  staff  in  both  developing  codes  and  running  problems .  The  machine 
used  in  the  beginning  was  the  Model  II  CPC  with  five  941  units .  The  limi¬ 
tations  of  this  equipment  forced  changes  in  methodology  which  later  proved 
applicable  to  codes  for  the  IM  701.  The  relative  power  of  the  final  701 
set-up  aroused  the  interest  of  IBM's  Data  Processing  Center  and,  with  the 
advent  of  the  704,  RAND  and  IBM  decided  to  develop  a  linear  programming 
system  for  this  machine  jointly.  Mr.  Hal  Judd  of  IBM  has  worked  with  Miss 
Cutler  and  the  writer  on  this  project  and  much  of  the  material  below  was 
originally  prepared  to  describe  the  704  system.  The  ideas  here  presented 
are  not  original  to  the  704  codes,  however,  auid  the  methods  are  not  tied 
to  any  one  machine,  except  as  noted.  Applications  to  other  machines  are 
discussed  briefly  in  Part  C. 
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The  growing  literature  on  the  subject  of  linear  programming  falls  roughly 
into  three  categories;*  (i)  applications  of  the  linear  programming  technique 
to  problems  of  the  rnildtary,  industry  and  business,  and  occasionally  to  prob¬ 
lems  simply  of  academic  interest;  (ii)  methods  for  solving  the  central  radthe- 
matical  problem  posed  by  the  resulting  models;  and  (iii)  computing  techniques 
for  implementing  these  methods.  Our  concern  is  mainly  with  the  last  although 
it  is  impossible  to  completely  separate  these  three  aspects.  A  concept  of 
analysis  has  value  only  if  it  shows  promise  of  being  applicable  to  real-lffe 
problems.  V/here  calculations  are  involved,  they  must -at  least  be  tractable 
in  the  abstract  and  feasible  within  current  computational  capabilities. 

It  is  not  Our  purpose  to  attempt  either  a  definition  of  linear  program¬ 
ming  and  its  uses  or  a  comparison  of  the  basically  different  methods  of  calcu¬ 
lation,  It  is  now  pretty  generally  agreed  that  the  simplex  method  of  Dr.  George 
Dantzig  provides  the  best  known  approach.  We  must  insist,  however,  that  it  is 
the  fundamental  theorems  of  this  method  which  are  of  importance  and  not  the 
particular  algorithm  and  tableau  which  were  set  down  in  the  original  exposi¬ 
tion.  The  calculations  can  be  carried  out  in  different  ways  and  adaptatibns 
and  extensions  are  es-sy  to  devise,  at  least  in  theory.  One  of  the  earlieSt 
and  most  widely  used  special  cases  -  the  transportation  problem  -  was  suffici¬ 
ently  important  to  warrant  development  of  efficient  computer  codes  for  this 
one  purpose,  but  it  is  difficult  to  find  other  sub-problems  as  well  definjed 
or  as  easily  simplified.  Consequently,  a  general  algorithm  has  been  applied 
in  nearly  all  other  cases  bvit  both  theory  and  experience  indicate  that  there 
is  an  upper  bound  well  below  500  to  the  number  of  restraints  that  can  be 

-((.  lieu  of  references  in  the  text,  a  catalogued  bibliography  is  givdn 

at  the  end  of  the  paper. 
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handled  efficiently  in  this  manner  with  foreseeable  machines.  Though  this 
figure  is  approximate,  it  is  probably  already  too  high  for  some  problems. 

This  has  been  the  cause  of  much  concern.  There  has  been  a  continuous  effort 
at  RAND  for  three  years  to  develop  general  codes  for  larger  and  larger  prob¬ 
lems.  Being  now  able  to  handle  some  250  restraint  equations  in  any  number 
of  variables  and  not  yet  having  satisfied  the  demands  of  many  realistic  prob¬ 
lems,  we  conclude  that  further  attempts  in  this  direction  will  be  fniitless. 

We  feel,  however,  that  the  present  codes  are  of  considerable  value  and,  more 
importantly,  that  their  evolution  has  shown  the  direction  which  future  develop 
ments  must  take.  In  this  regard,  perhaps  the  most  important  result  is  the 
mutual  realization  by  those  who  formulate  problems  and  methods  and  those  who 
devise  computer  programs  that,  for  better  or  worse,  they  are  joined  together. 

Our  aim  is  to  discuss  the  growth  of  these  codes,  from  the  standpoint  of 
both  mathematical  and  coding  techniques,  and  then  to  give  a  brief  resume  of 
current  and  proposed  developments  in  this  field  at  RAND.  The  presentation 
he  in  three  main  parts  which  will  bear,  more  or  less  respectively,  on 
the  following  three  theses. 

A.  In  order  to  handle  efficiently  a  large  number  of  problems  of  one 

type,  codes  must  be  designed  with  an  eye  to  the  two  aides  of  com-  ' 
putational  problems:  (1)  the  method  must  be  adapted  to  the  char¬ 
acteristics  of  the  machine,  and  (2)  provision  must  be  made  for 
handling  easily  any  reasonable  special,  requirements,  for  making 
a  variety  of  decisions  automatically  and  for  easy  operation. 

B.  In  order  to  attain  the  required  efficiency  and,  flexibility,  the 
codes  must  themselves  be  organized  in  a  manner  analagous  to  the 
method  and  not  based  on  standard  routines  or  universal  methods. 

C.  For  solving  very  large  problems  at  reasonable  cost,  it  will  be 
necessary  to  develop  special  techniques  for  special  classes. 
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PART  A  -  THE  REVISED  SIMPLM  METHOD  WITH  SPECIAL  FEATURES 

A  general  familiarity  with  the  simplex  method  will  be  aisswned.  For 
example,  no  proofs  will  be  given  for  the  standard  theorems  concerning  basic 
solutions  or  the  simplex  criterion.  However,  the  algorithm  using  the  pro¬ 
duct  form  of  the  inverse  and  certain  other  features  will  be  developed  in 
detail.  The  Part  is  divided  into  twelve  sections.  Sections  X  and  XII  con¬ 
tain  descriptions  of  devices  not  previously  written  up  in  detail  but  they 
have  been  in  use  for  some  time  and  their  usefulness  is  established. 

I  -  EARLY  CPC  SET-UPS  AT  RAND 

Since  it  was  felt  that  the  transformation  of  the  entire  tableau  in  the 
original  simplex  method  caused  unnecessary  work  and  a  waste  of  storage  capacity, 
especially  in  problems  with  many  more  variables  than  restraints,  Dantzig 
proposed  that  the  CPC  be  set  up  for  the  revised  simplex  method.  Here  only 
the  inverse  of  the  basis  is  maintained,  together  with  the  transformed  right 
hand  side  and  a  pricing  vector  which  form  an  extra  column  and  row.  Since 
examples  of  "cycling"  due  to  degeneracy  are  extremely  rare  and  artificial, 
no  provision  was  to  be  made  for  rigorously  resolving  this  matter.  Previous 
experience  with  hand  calculations  and  on  the  SEAC  had  indicated  that  such  a 
decision  caused  no  practical  difficulty. 

It  was  realized  that  the  CPC  was  an  inadequate  machine  for  linear  pro¬ 
gramming  but  it  was  thought  that  it  would  provide  a  good  training  and  test¬ 
ing  ground  during  the  interim  until  large  machines  were  available.  This 
was  borne  out  by  subsequent  experience. 

After  some  months  of  study,  planning  and  board  wiring,  a  set-up  was  in 
operation  which  would  handle  2?  restraints  in  70  variables.  Operation  was- 
fairly  automatic,  it  being  necessary  to  feed  through  two  decks  of  cards  for 
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each  iteration.  But  it  was  painfully  slow,  due  mainly  to  the  punching  of 
a  new  inverse  deck  on  each  cycle.  V/atching  this  operation  one  time  was 
enough  to  convince  anyone  that  the  whole  set-up  should  be  junked. 

Although  most  elements  of  the  inverse  were  zer^  during  the  first  few 
iterations,  the  nxamber  of  non-zero  elements  increased  on  succeeding  cycles. 

It  was  awkward  to  delete  zeros  since  the  necessary  bookkeeping  really  demand¬ 
ed  a  stored  program  machine.  However,  even’ on  a  large  machine,  a  collapsed 
inverse  matrix  requires  excessive  bookkeeping  and  searching  if  both  the 
inverse  and  its  transpose  are  needed,  that  is,  if  sometimes  one  wishes  combina¬ 
tions  of  rows  and  sometimes  combinations  of  columns.  Although  this  could 
be  avoided  in  the  straight  simplex  method,  it  is  the  kind  of  flexibility  that 
is  needed  for  variations. 

Indeed,  another  shortcoming  of  the  set-up  was  the  very  fact  that  the  CPC 
had  been  made  into  a  'black  box*'  which  rendered  it  inflexible  and  hence  im¬ 
practical  in  use.  Model-makers,  like  women,  reserve  the  right  to  change  their 
minds . 

Dantzig  recalled  an  old  suggestion  of  Alex  Orden  on  the  product  form  of 
inverse  and  upon  investigation  this  approach  turned  out  to  resolve  most  of  the 
difficulties  then  current.  A  new  CPC  set-up  was  designed  which  coTild  handle 
40  restraints  in  99  variables  (or  more  with  a  little  care)  and  was  used  suc¬ 
cessfully  on  several  problems.  Among  these  was  a  series  of  models  in  a  study 
of  petroleum  refinery  scheduling  being  conducted  by  Dr.  Alan  S.  Manne  and 
economic  processing  models  in  general  being  conducted  by  Dr.  Harry  Markowitz. 
More  work  was  imposed  on  the  machine  operator  but  considerable  flexibility 
was  attained.  Based  on  work  of  Manne,  Dantzig,  Markowitz  and  the  writer, 
techniques  were  arranged  to  vary  the  right  hand  side,  the  optimizing  form  and 
single  elements  of  the  matrix.  These  are  described  in  later  sections. 
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II  -  THE  RAND  SIMPLEX  CODES  FOR  THE  701 

In  mid-1953  work  was  begun  on  a  701  code  to  employ  the  general  features 
of  the  second  CPC  set-up.  As  many  are  well  aware,  preparation  for  the  701 
was  a  big  step  from  the  CPC,  regardless  of  the  problems  worked  on.  Many  mis¬ 
takes  were  made  and  later  corrected  as  sophistication  was  gained.  However, 
the  writer  fell  into  one  fallacy  which  is  hard  to  explain,  let  alone  justify. 

The  CPC  set-ups  had  used  8-diglt  floating  decimal  arithmetic.  As  many 
as  80  iterations  on  35-order  systems  had  been  performed  before  accumulated 
round-off  error  became  troublesome.  The  initial  goal  for  the  701  was  100 
restraint  equations  and  high  operating  speed.  In  trying  for  speed,  it  was 
decided  to  use  single-precision,  fixed-point  arithmetic.  The  code  was  tested 
on  a  24-order  system  consisting  mainly  of  the  identity  matrix.  But  when  a 
dense  7-order  system  was  tried,  the  code  stopped  after  a  few  iterations  with 
an  error  indication.  After  some  hand  computations,  it  was  realized  that  the 
error  stop  was  legitimate,  it  was  actually  due  to  round-off  1  The  code  was 
revised  in  an  attempt  to  squeeze  out  more  significant  places,  but  to  no  avail. 
The  basic  transformations  were  simply  too  inaccurate. 

Finally  realizing  the  significance  of  the  CPC  experience,  we  rewrote  the 
transformation  sections  of  the  code  to  perform  double-precision,  floating¬ 
point  arithmetic,  i.e.  carrying  nearly  16  decimal  digits  (52  binary  bits)  pl\is 
a  binary  exponent.  This  permitted  100-order  systems  to  be  solved  and,  in  fact, 
250  equations  are  now  handled  with  essentially  the  same  arithmetic  (54  binary 
bits  plus  exponent.)  It  has  been  possible  to  keep  running  time  down  on  some 
operations  by  using  fixed-point  arithmetic,  partly  double  and  partly  single 
precision.  We  feel  that  the  added  complexity  is  a  small  price  to  pay  for  the 
results.  The  inputs  are  still  required  to  fall  into  a  restricted  fixed-point 
range.  This  causes  no  difficulty  in  linear  systems  usually  of  low  precision 


and  ia  an  .ffeotiva  way  to  obtain  reasonably  good  scaling  trom  a  variety  of 
users  of  the  system. 

Other  more  excusable  flaws  began  to  appear.  The  initial  flexibility 
attained  on  the  CPC  was  discouragingly  difficult  on  the  701.  Improvements 
and  changes  for  easier  operation  or  special  requirements  involved  prodigous 
tasks  of  code-checking.  Further,  some  believed  that  a  high  speed  computer 
code  should  read  original  data  cards,  solve  the  problem,  print  the  answers 
and  stop.  We  tried  this  mode  of  operation  and  found  it  impractical  for  large 
computations  with  large  amounts  of  input.  For  reasons  well  known  at  any  com¬ 
puting  center,  the  probability  is  practically  zero  that  a  large  Job  will  go 
to  completion  with  no  error  and  in  one  run.  There  is  the  added  uncertainty 
as  to  whether  the  right  problem  got  into  the  macldne  correctly  in  the  first 

place.  We  have  had  the  experience  of  getting  correct  answers  to  the  wrong, 
probl em . 

Consequently ,  the  whole  program  was  completely  re-organized  and  divided 
into  separate  parts.  The  initial  conversion,  editing  and  packing  of  data  is 
now  done  by  a  separate  data  assembly  code  which  produces  binary  cards  and 
magnetic  tape.  Errors  are  expected  and  provisions  made  beforehand  for  recov¬ 
ery.  The  main  code  is  organized  into  replaceable  or  interchangeable  regions 
which  correspond  to  the  major  divisions  of  the  method.  The  present  programs, ’ 

which  are  the  result  of  two  subsequent  revisions,  are  discussed  further  in 
Part  B. 

Ill  -  ADVANTAGES  OF  THE  PRODUCT  FORM  OF  INVERSE 

The  advantages  claimed  for  the  revised  simplex  method  with  the  product 
form  of  inverse  will  now  be  summarized  and  the  algorithm  then  developed  in 
detail.  It  is  necessary  to  keep  in  mind  the  following  facts  concerning  the 
matrix  of  coefficients  of  a  typical  linear  programming  problan  and  the  re- 
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quirements  of  the  simplex  method o 

(a)  The  matrix  is  usually  very  sparse,  the  percentage  of  non-zero 
elements  being  typically  5  to  15  per  cent  of  the  total. 

(b)  The  given  data  is  almost  invariably  of  low  precision  and  easily 
scaled. 

(c)  Some  problems  have  a  large  ratio  of  number  of  variables  to  niimber 

of  restraints,  between  5  to  1  and  10  to  1  being  fairly  common.  (One 
problem  run  on  the  701  had  a  ratio  of  about  30  to  1.) 

(d)  The  rules  of  elimination  in  effecting  a  change  of  basis  (where  two 
successive  bases  differ  only  by  one  colximn)  are  the  same  no  matter 
what  particular  algorithm  is  used.  Although  such  a  transfonnation 
is  easily  represented,  its  total  effect  on  the  numerical  represen¬ 
tation  of  a  set  of  vectors,  i.e.  a  matrix,  may  be  extensive. 

V^ith  these  observations  in  mind,  the  advantages  of  the  present  method 
can  be  stated  as  follows. 

(1)  Since  the,  original  data  matrix  is  not  transformed  from  iteration 

to  iteration,  it  is  clear  from  (a)  and  (b)  that  an  elaborate  organ¬ 
ization  of  the  data  can  be  performed  at  the  outset,  once  and  for 
all,  so  that  it  is  in  the  most  convenient  and  compact  form  for  use 
throughout  the  problem.  Since  the  original  matrix  is  referred  to 
only  once  during  the  cycle,  it  can  be  stored  in  auxiliary  storage 
(say  on  a  tape)  out  ox  the  way,  and  its  compact  form  minimizes 
transfer  time  as  well  as  computing  time. 

(2)  It  is  clear  from  (c)  that  it  would  be  very  expensive  in  some  prob¬ 
lems  to  transform  the  whole  matrix  on  each  cycle.  If  there  are  n 
variables  and  m  constraints  and  if  n  is,  say,  only  5ni,  it  is  very 
unusual  for  the  process  to  take  5m  iterations.  In  other  words. 
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many  vectors  (colxunns)  are  never  selected  at  all  to  come  into  the 
basis  and  hence  there  is  no  use  doing  anything  to  them. 

(3)  As  remarked  above,  the  modification  and  re-recording  of  the  m^  ele¬ 
ments  of  a  full  inverse  is  awkward  and  time-consuming,  especially 
if  zeros  are  deletedo  On  the  other  hand,  the  produce  form  requires 
the  recording  of  only  one  additional  column  (plus  an  index)  on  each 
iteration.  These  columns  will  ailmost  surely  contain  many  zeros  and 
may  be  condensed  since  they  are  always  applied  column-wise  and  re¬ 
cursively  to  a  single  Sector. 

(4)  The  product  form  is  extremely  amenable  to  modifications  of  the  method 
since  the  inverse  and  its  transpose  are  equally  easy  to  apply. 

One  apparent  disadvantage  of  the  product  fom  is  that,  while  only  one  new 
coluxnn  is  recorded  each  cycle,  it  is  necessary  to  read  T  columns  to  apply  this 
form  where  T  is  the  number  of  iterations  already  performed o  Thus,  after  ra 
iterations,  where  m  is  the  number  of  restraints,  one  must  read  more  information 
than  in  reading  a  full  inverse  (not  taking  condensation  into  account.)  This 
is  still  profitable  for  something  over  2m  iterations  as  can  be  seen  as  follows. 

With  an  explicit  inverse,  we  must  read  it  twice  and  write  it  once 
each  iteration,  giving  3niT  columns  handled  on  T  iterations., 

With  the  product  form,  we  must  read  t-1  columns  twice  on  iteration 
t  and  write  one  new  column,  giving 

T 

=  T(T-1)  +  T  “  T^  columns  handled  on  T  iterations o 

t*l 

2 

Setting  T  =  3niT  gives  T  »  3mo  However  something  must  be  allowed 
for  the  fact  that  less  arithmetic  accompanies  the  writing  than  the  reading o 
This  disadvantage  is  more  apparent  than  real,  however.  After,  say,  2m  itera¬ 
tions,  roxmd-off  error  will  begin  to  be  noticeable  on  large  problems  no  matter 
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which  form  is  used.  Many  problems  are  solved  in  2m  iterations  or  less  but,  ' 
if  it  takes  more,  one  can  re-invert  the  basis  matrix,  at  any  time,  producing 
not  more  than  m  colximns  of  information.  The  time  for  this  inversion’  is  much 
less  and  the  accuracy  of  the  resulting  transformations  is  as  good  or  better 
than  after  the  same  number  of  full  simplex  cycles.  (The  order  of  elimination 
in  inverting  is  designed  to  maintain  accuracy.)  A  special  code  is  provided 
for  this  purpose.  It  is  useful  for  solving  any  system  of  linear  eq\iatiops, 
especially  v/here  several  right  hand  sides  are  used  with  a  sparse  matrix. 

It  is  helpful  with  this  form  of  inversion,  to  consider  the  matrix  of 
coefficients  of  the  restraint  equations  as  a  collection,  or  set,  of  column 
vectors,  ignoring  the  fortuitous  ordering  given  to  the  variables  in  the  for¬ 
mulation  of  the  problem.  Whatever  scrambling  of  columns  that  may  occur  in 
the  process  is  recorded  in  a  list  of  basis  headings  which  accompany  and  iden¬ 
tify  the  basic  variables  of  a  particular  solution.  This  point  of  view  will 
be  adopted  throughout  the  present  discussion. 

IV  -  NOTATION 

In  discussions  of  matrices,  their  transposes,  columns,  rows,  elements, 
etc.,  there  has  arisen  in  the  literature  a  hodge-podge  of  notation  -  .upper- 
and  lower-case,  bars,  stars,  primes,  underscoring,  varied  type  faces,  etc. 
Hence  we  introduce  the  following  simplified  and  unified  notation  aimed  at 
reducing  the  number  of  symbols  and  type  forms  required  to  talk  about  matters 
of  inherent  simplicity. 

Real  nvimbers  will  be  represented  by  lower-case,  light  face  Latin  or  Greek 
letters.  However,  these  letters  will  be  indexed  as  necessaiy  and  an  indexed 
quantity  may  mean  one  particiilar  element,  a  row,  a  column  or  a  matrix.  The 
context  will  always  make  clear  which  is  meant.  In  exceptional  cases,  specific 
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remarks  will  be  made  on  the  side. 

The  distinction  between  Latin  and  Greek  letters  will  appear  latero  In 
either  case,  a  row  vector  will  be  denoted  by  a  general  subscript,  where  the 
range  of  values  which  this  index  takes  on  is  defined  when  it  is,^  first  men¬ 
tioned,  An  element  of  the  row  vector  will  be  denoted  by  a  specific  subscript. 

c  is  a  row  Mrector,  c-  is  the  element  of  index  3« 

A  column  vector  will  be  denoted  by  a  general  superscript,  one  of  its  elements 


by  a  specific  superscript. 

-la  a  n 


is  a  colxmin  vector. 


is  the  element  of  index  2. 


A  matrix  will  be  denoted  by  both  a  superscript  and  a  subscript, 
a^  is  a  matrix 


is  the  row  consisting  of  t  he  element  of  index  5  in  each  column. 


a^  is  the  column  consisting  of  the  element  of  index  4  in  each  row. 
4 

a^  is  ihe  element  in  row  6  and  column  2. 

It  is  important  to  note  that  no  sifnificance  is  attached  to  the  letter 
used  for  an  index.  Thus,  a^  and  refer  to  the  same  quantities.  However, 
a^  is  the  transpose  of  a^.  The  transpose  of  a  matrix  must  be  denoted  by  defin¬ 
ing  a  new  letter,  e.g,,  b^  -  a^  defines  b^  as  the  transpose  of  a^. 

We  will  also  adopt  a  modified  form  of  the  summation  convention,  as  follows. 
Siitmnntion  Convention;  Whenever  the  same  letter  appears  twice  in  a  term, 
as  the  index  of  a  row  and  then  following  as  the  index  of  a  column, -this 
letter  is  a  dummy  index  and  summation  is  understood  over  the  (previously 
defined)  range  of  these  indices. 


eous 


For  given  coefficients  and  constants  b^,  then,  a  system  of  m  simultan- 
linear  equations  in  n  unknowns  x?  would  be  written  simply 


(i  =  1,  o..,  m;  j  ”  1> 
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n). 


A  matrix  may  be  multiplied  on  the  left  by  a  row  vector  to  produce  a  row  vector: 

If  a  matrix  is  multiplied  on  the  left  by  a  row  and  on  the  right  by  a  column, 
the  resxilt  is  a  scalar: 

c.  a  j  x"^  *2  . 

Note  that  a^  b^^  does  not  denote  summation.  If  it  means  anything,  it  simply 

denotes  a  set  of  products.  The  quantity 

i  i  , 
c  .  =  a  b  . 

'  '  i 

is  a  matrix,  the  so-called  outer  product  of  a  and  b^. 

Uppercase  letters  (except  T)  v/ill  be  used  for  sets.  The  use  of  letters 
as  indices  is  as  follows; 


h,  i,  j,  foi'  general  indices, 

m,  n,  for  limits, 

p,  q,  r,  3  for  specific  indices. 

When  it  is  desired  to  index  a  quantity  over  time,  i.e.  iterations,  an  index 
(denoted  in  general  by  t)  enclosed  in  parentheses  will  be  used,  e.g. 
is  produced  from  after  t-1  transformations.  The  capital  letter  T  will 

s(ij 

be  used  as  the  current  limit  for  t,  that  is,  T  is  the  current  iteration  n’om- 
ber  and  t*l,  ...,  T. 


V  -  FORM  OF  THE  PROBLEM  AMD  INITIAL  SET-UPS 

It  is  mandatory  to  start  with  the  identity  matrix  as  the  initial 
basis  in  all  circumstances,  we  define  the  standard  linear  programming  problem 
in  two  forms.  Certain  deviations  from  these  standard  forms  are  discussed  below. 


*  Note  the  similarity  to  the  notation  Ax  - 
and  b.  Our  summation  convention  is  simply  the 
cation  combined  with  indicial  notation. 


b  for  matrix  A  and  columns  x 
usual  rule  for  matrix  •mxlLtipli- 


' '  ^ 
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Standard  Form  1 

Given;  A  row  of  coefficients  of  a  linear  fornj  to  be  optimized 
a°  (j  -  1,  2,  i) 

A  matrix  of  restraint  coefficients 

Ej  (i  =  1>  2,  ...»  m;  j  “  1»  2,  ...»  l') 

A  column  of  constants 

(i  “  1»  2,  .. .»  m) 

To  find:  A  column  of  values  for  x*^  >  0  (j  =  1,  2,  ..o»  2)  such  that 
the  variable  x°  is  maximized  subject  to 

(5.11)  +  a?  =  0 


(5.12)  a^  x-^  <  b^  . 

Note  that  (5.11)  is  perfectly  general  since  the  siim  a^  x^  inay  be  minimized  or 
maximized  merely  bj’"  changing  the  signs  of  the  a^  .  To  convert  (5.12)  to  equal¬ 
ities »  we  define  the  variables  x'^  ^  >  0  (i  =  1»  ...»  m).  Then  (5 <.12)  becomes 


(5.12')  aj  xJ  +  =  b^  . 

If  c)^  (i’-0,l»  ...»  m;  h-0,l»  ...»  m)  is  the  (m-fl)-order  identity  matrix 

(essentially  the  Kronecker  delta) »  then  we  can  define 
a^  =  Sq  and  b*^  “  0 


thus  defining  j  =  0  and  j  =■  2+l»  ...»  2+m  «  n»  and  replace  both  (5oll)  and 


(5cl2')  by 

(5.13)  a^  x^  =>  b^  (i  =  0»  1,  ...»  mj  j  “  0»  1»  ...»  n). 

t3  , 

If  the  restraint  equations  can  be  put  in  the  form  (5.12' )»  then  some  of  the  b 
are  permitted  to  be  nQgative»  if  necessary.  (See  section  X  below.)  However 
the  number  of  b^  <  0  should  be  small  to  prevent  an  excessive  number  of  itera¬ 
tions.  If  the  restraints  cannot  be  specified  exactly  in  the  form  (5.12) »  then 
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the  problem  should  be  put  in  Standard  Form  2. 

Standard  Form  2 

Given;  A  row  of  coefficients  of  a  linear  form  to  be  optimized 

(j  «.  1,  2,  n) 

A  matrix  of  restraint  coefficients 

(1  “  2,  3,  mj  j  =  1,  2,  n) 

A  column  of  constants 

b^  >  0  (i  =  2,  3,  ...»  m) 

To  find:  A  column  of  values  for  >  0  (j  =  1,  n)  such  that  the 

variable  is  maiximizea  subject  to 

(5.21)  x^  +  a^  x*^'  “  0 

aJxJ  -  (1  .  2,  3,  «). 

We  now  add  to  the  system,  artificially,  the  identity  matrix  (except  5^ 
which  is  always  incorporated)  and  auxiliary  variables  (k  -  1,  ...,  m) 

and  construct  an  auxiliary  form  in  which  the  variable  is  to  be  maximized 

first,  i.e.  before  maximizing  x°.  This  process  is  called  Phase  I. 

The  purpose  of  Phase  I  is  to  eliminate  the  artificially  added  columns 
from  “the  system  or  at  least  to  make  sure  that  the  corresponding  variables 
become  zero.*^  Hence  we  put  a  weight  of  unity  on  each  one  in  the  row  of  index 
1,  which  will  become  the  auxiliary  form.  Let  a^^^  =  and  for  k  «  2,  3,  . . . ,m 

^n+k  "  ^1  ^  ^k  *  Variables  x  are  to  be  non-negative  except  for 

x’^^^  which  is  defined  by 

(5.23)  x"-^^  ^  -  0  . 

k-2 

At  least  one  artificial  column ,  usually  the  column,  must  remain  in  the 
solution  (basis)  at  zero  level  even  while  maximizing  the  variable  x^. 
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the  row  1  equation  takes  the  form 

(5.27)  a^xJ  +  (j  =  0,  1,  n). 

i  i  •  4^ 

Furthermore  the  columns  a,';,  a  ,,  a  ,,,  nOw  form  . 

0  n+i'  n+2  ’  n+ra  “h 


The  data  assembly  program  computes  a^  and  automatically.  If  some, 
but  not  all,  of  the  columns  of  occur  in  (j<n),  then  these  columns  should 
be  indicated  as  being  in  the  initial  basis o  The  corresponding  rows  will 


then  be  omitted  in  the  sums  (5*26) o  This  is  equivalent  to  adding  artificial 
columns  only  for  those  columns  of  which  are  missing. 

The  vectors  a^,  never  entered  with  the  data.  They  are  implicit 

in  the  code  and,  once  eliminated,  cease  to  exist.  The  basis  headings,  which 
are  the  J-indices  of  the  columns  a^  in  position  h  =  0,  1,  ..o,  m  of  the 
basis,  are  left  zero  for  all  artificial  vectors  since  the  position  h  identi¬ 
fies  them  sufficiently.  Legitimate  columns  of  on  the  other  hand,  may  go 
out  of  the  basis  and  come  back  in  out  of  position.  Hence  they  must  have  names. 


as  in  Standard  Form  1<, 


The  colxjmn  indices  j  =  1,  ...,  n  are  used  above  for  expository  purposes 

only.  Any  n  distinct  symbols  can  be  used  since  they  are  only  names,  no  parti¬ 
cular  ordering  is  implied.  The  same  is  not  true,  of  course,  of  the  row  indices. 


Alternate  Methods  of  Starting  a  Problem 

Experience  has  led  to  the  incorporation  of  another  device  for  obtaining 


initial  solutions.  Sometimes  the  formulator  of  the  problem  knows  a  feasible 


basis  other  than  the  identity.  Provision  is  made  for  introducing,  arbitrarily, 
any  number  of  column  vectors  into  the  basis  at  the  outset,  with  the  machine 
making  the  decision  as  to  which  column  of  the  basis  each  should  occupy.  If 
the  formulator  has  misjudged  and  the  specified  columns  produce  a  8ingul.ar 
matrix,  the  machine  prints  out  (and  saves)  the  pertinent  information  and 
stops.  If  the  resulting  matrix  is  non-singiolar  but  the  solution  is  infeas- 
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ible,  then  the  composite  algorithm  automatically  cuts  in  and  works  toward 
feasibility  in  succeeding  transformations.  (See  section  X  below.) 

One  other  provision  of  a  less  absolute  nature  has  proved  useful.  Occa¬ 
sionally,  a  model  is  a  re-work  of  an  older  problem  so  that  something  is  known 
about  its  behavior.  At  other  times  the  formulator  has  certain  insights  which 
he  would  like  to  exploit  without  committing  himself  to  absolute  statements 
regarding  feasibility  or  singularity.  It  is  possible  to  arrange  the  columns 
of  aj  in  order  of  decreasing  likelihood  of  use,  that  is,  with  the  most  like¬ 
ly  candidates  for  entrj^  into  a  feasible  or  optimal  basis  occuring  first. 

These  may  be  separated  from  the  others  by  a  special  mark  (called  a  "  curtain"  ) 
which  is  equivalent  to  the  following  instruction:  "  If  any  candidate  for  entry 
into  the  basis  is  available  ahead  of  this  mark,  use  itj  otherwise  proceed  to 
the  others."  Several  marks  may  be  used.  They  have  often  reduced  the  number 
of  iterations  required  but  when  used  injudiciously,  i.e.  on  a  mere  hunch,  may 
have  the  opposite  effect. 

There  are,  of  course,  other  devices  which  require  no  special  provision^ 
except  perhaps  for  loading.  For  example,  certain  activities  (columns)  may 
be  withheld  from  the  machine  until  optimality  is  obtained  with  the  others. 

This  is  another  advantage  of  the  revised  metJiod  and  additional  data  may  be 
added  or  subjected  to  the  optimality  test  at  any  time,  simply  by  re-loading. 


71  -  THE  BASIS  AND  ITS  INVERSE 


be  denoted  by 

o«.»  nij  h  «  0,  1,  0..,  mj  °  ^0  ^ 

where  it  is  always  assumed  that  pj  ==  aj  =  .  The  matrix  pj  is  called  the 


basis  and  it  changes  from  iteration  to  itBration„  One  always  begins  the  com- 


putations  with  the  initial  basis 
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. 

h 


This  imples  that  a.  (perhaps  augmented  artificially  as  in  Standard  Form  2) 

•J 

i.  i  i 

•  contains  .  Since  t  he  columns  of  Sj.j  ®ay  occur  in  any  order  in  a^  ,  we  use 

second-order  subscripts  on  the  j*s  to  denote  these  column  names,  i.e. 


(6.01) 


Initially. 


if  T  ^  i 

After  T  cycles,  some  of  the  colximns  of  ^  will  differ  from  .  The 
coltamns  have  been  drawn  from  but  if  we  write 
(6.02) 

the  in  (6.02)  are  not  the  same  as  in  (6.01).  Hence  it  is  necessary  to 

tag  thie  with  T  and  keep  the  definition  of  up  to-date.  These  indices 

are  called  basis  headings  and  are  referred  to  frequently. 

i(T) 

We  will  not  be  so  much  concerned  with  the  basis  '  on  iteration  T  as 

i(T) 

with  its  inverse  whiich  will  be  denoted  by  ,  that  is 

(L  i(T)  k(T)  i(T)  k(T)  ci 

The  matrix  never  exists  explicitly  but  is  carried  as  a  product  of  ele- 

h 

mentary  matrices  which  are  stored  in  a  condensed  form  consisting  of  at  most 
one  column  plus  an  index  and  identification.  These  columns  are  t  heraselves 
condensed  with  only  non-zero  elements  recorded.  They  are  called  t rang f omat i on 
vectors  (sometimes  elimination  'Equations"  ,  see  (6.5))  and  denoted  by 
Since  the  index  r  is  itself  a  fxinction  of  t,  it  should  be  written  r(t).  However 
this  will  not  usually  be  done  since  it  sometimes  appears  as  a  superscript  in  an 
array  already  indexed  (t-1)  and  confusion  would  result. 

As  already  indicated,  Greek  letters  will  be  used  for  a  basis  and  its 
inverse,  and  also  for  all  (m+1) -order  vectors  expressed  in  terms  of  the  basis. 
Likewise  certain  ratios  and  functions  involved  in  decisions  concerning  the 
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basis  will  be  denoted  by  Greek  letters. 

rati 
i(0) 


An  explanation  of  the  generation  and  use  of  the  vectors  is  in  order. 


i(t) 

1  '  ‘  V€ 


Since  we  started  with 


on  hand  to  start  iteration  T.  Suppose  that  column  a 


for  the  first  iteration,  we  will  have 
i 


s(T) 


i  ifT— 1’)  ifTl 

a^  to  replace  column  h  =  r(T)  in  producing  a  new  basis  ■■  ^  ‘ 


h 

has  been  chosen  from 

p^'  • .  Through¬ 


out  the  following  discussion,  the  specific  index  r  is  to  be  understood  as  r(T) o 
Superscripts  in  parentheses  will  always  refer  to  the  entire  array  to  which  the 
element  belongs  and  are  not  to  be  understood  as  modifying  the  open  superscript. 
Let  satisfy  the  equation 


(6.2) 

i 


i(T-l)  h 


s(T)  s(T)  * 


Clearly  can  be  obtained  by 


(6.3) 


i(T-l)  h 
"h  ^s(T) 


Now 


“s(T)  • 

i(T)  ^k  ^  j_i?-(T“l)  foj.  h  ^ 


(6.4)  pj^-  •  0^  -  pj^- 

since  only  column  h  =  r(T)  changes.  However,  from  (6.2)  we  can  solve  for  the 


exceptional  column. 

(6.5) 


Vt) 


i  ^  i(T-l)  h 
^(T)  "  ^"h  ^ 


(h  r  ) 


where 

(6.6) 

Now  defining 

(6.7) 


-a 


sill 


3(T) 


(  /  0  by  choice)  . 


r(T) 


“s(T) 


i(T) 


(i  r) 


i(T) 

% 


(h  /  r) 


we  can  replace  (6.4)  and  (6.5)  by 

(6.8) 


i(T)  k(T)  „  i(T-l) 


i  1(T)  ifT^ 

since  becomes  column  p^'  ' .  Clearly  only  the  column  ' 


and  the  index 


r  ■  r(T)  need  be  recorded.  Now  multiplying  both  members  of  (6,8)  on  the  left 


by  and  on  the  right  by 


we  obtain 
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.  (6.9)  .  .J(T)  . 

Equation  (6.9)  is  the  heart  of  the  product  form  of  inverse.  Applying  it  for 
t  “  1,  2,  T  and  using  h^  to  indicate  dummy  indices,  we  obtain 


(6.10)  21’^)  . 


^*^T-1  ^^T-2 


^h 


t-1 


^h 


Hence  an  equation  like  (6.3)  implies  the  recursive  generation  of  its  right 
member  by  using  the  form  of  given  by  (6.10)  „  This  is  easily  done  as 


follows,  using  (6.3)  for  an  example. 

^s(T)  “ 


Let 

form 


“s(l) 


i(l)  -h  -i 

^h  “s(l)  “  %(2) 


(6.11) 


i(t-l)  -h  -i 

'^sCt-l)  “  “ 


3(t) 


i(T-l)  -h  -i  i 

’^h  “s(T-l)  ”  “s(T)  °  “3(T) 

“s(t-i;  “3(t-l)  (i  /  r  -  r(t-iy  ) 

®3(t-l)  (^  “  r(t-l)  ) 

suffices  for  multiplying  Tt^('^^  by  a  row  vector  on  the 
left.  Suppose  it  ie  necessary  to  compute 
(6.13)  ,J(T-1)  .  ,(l)  _ 

We  use  the  transformation  vectors  in  reverse  order  to  that  of  (6.11)  ; 

Let  f.  =" 

1  n 

-d)  h(T-l)  _  -(2) 


It  is  easy  to  see  that 
-i 

a  \  " 
s(t) 

(6.12) 

“s(t)  “ 

An  equally  simple  rule 


(6.14) 


form 
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(6.14  cont.) 


-(T-1)  h(l)  -(T)  (T) 


Again  it  is  easy  to  see  that 


(6.15) 


(1  /  r(T-t)  ) 


i(T-t) 


VII  -  THE  PRICING  OPERATION  rchoosin/r  index  s') 

The  row  vector  in  (6.13)  is  called  the  pricing  vector.  In  the  simp- 

lest  problem  (Phase  II,  maximize  x^,  no  infeasibilities),  .  However, 

in  a  typical  Phase  I,  =■  ,  and  in  general  f j,  =  1  for  several  i. 

The  pricing  vector  is  applied  to  a^  , 

(7.1)  »  a(T) 


(The  dj  are  what  are  called  (+)  (c^  -  z J  in  the  original  simplex  method.) 
To  choose  to  bring  into  the  basis,  take 


(7.2)  d^'^^  *  min  d^'^^  <  0  . 

®  J 


(If  the  minimum  is  not  unique,  the  first  such  index  s  is  retained.) 

If  all  dj  >  0,  the  phase  is  complete.  In  a  Phase  I,  this  is  Terminal  Ij 
in  a  Phase  II  it  is  called  Terminal  2  and  is  the  point  at  which  one  usually 
expects  to  arrive  and  quit,  i.e.  the  optimal  solution  is  attained. 

VIII  -  CHOOSING  THE  BASIS  COLUMN  TO  BE  REPIAHKn  (Index  r) 

Let  the  current  basis  be  the  current  optimizing  form  be  row  p 

and  the  current  solution  vector  be  i.e. 


(8.1) 


pi(T-l)  .  i(T-l)  ^h  ^ 


We  will  assiome  the  solution  is  feasible,  that  is’ 


(8.2) 


pi(T-l)  ^  Q  for  i  >  q 
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where  row  q  is  the  first  actual  re'straint  equation.  Normally  q  =»  p+1 . 

The  usual  criterion  for  choosing  r(T)  is  to  choose  follows, 

Let  A  be  the  set  of  indices  i  >  q  for  which  >  0*  Then 


(8.3) 


r(T)  ilA 


,i(T-l) 


“s(T) 


>  0  . 


i ( T-1 ) 

The  problem  of  degeneracy,  i.e.  multiple  values  of  i£  A  for  which  0  ^  =  0,  is 

disregarded  except  for  the  following  rule  which  has  proved  effective  for 
reducing  round-off  error.  (It  incidentally  breaks  the  tie  in  Hoffman's  examples 
of  "cycling."  )  If  *  0  and  r(T)  is  ambiguous,  choose  r(T)  so  that 

largest  possible  (positive)  value.  In  case  of  further  ambiguity, 
take  the  smallest  such  index. 

We  will  modify  (8.3)  slightly.  Let  R  be  the  set  of  indices  i  >  q  for 
which  and  have  the  same  sign  with  /  0*  Then  let 


(8.4) 


r(T)  if.  R 


“s(T) 


>  0  . 


Note  that  (8.4)  gives  the  same  result  as  (8.3)  as  long  as  (8.2)  holds.  Degen¬ 


eracy  is  resolved  in  the  same  way,  that  is,  if  9 


0,  take  oLfrni  >  0  and  max. 


If  no  0  ,  V  can  be  chosen,  that  is,  all  ratios  are  non-positive  and  zero 
r(T) 

ratios  have  negative  denominators,  then  has  no  finite  maximum.  A  class  of 
solutions  exist  as  follows: 


(8.5)  -  9 


with  the  value 

(8.6) 


,pCT-r)  .  9  „P 


+0)  as  0 


This  is  Terminal  3.  It  cannot  happen  in  a  Phase  I  since  clearly  zero  is  an 


upper  bound  for  the  variable  x 


£ 
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Clearly,  <  0  and,  if  ••  0,  then  all  “0,  When  (and  if)  this 

condition  is  attained,  the  artificial  variables  are  maintained  at  zero,  includ¬ 
ing  by  considering  (5.23)  as  a  restraint  while  maximizing  x^  (Phase  II.) 

If  x*^^^  cannot  be  driven  to  zero,  there  is  no  solution  to  the  given  problem. 
This  condition  is  called  Terminal  1. 

The  problem  at  this  point  can  be  displayed  in  the  follo'Aang  augmented 
form  where  the  are  shown  above  the  matrix  of  detached  coefficients. 
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IX  -  SUMMARY  OF  CYCLIC  OPERATIONS  (One  Iteration) 

(9.0)  Test  for  end  of  a  Phase  I  or  for  arbitrary  halt  (external  switch.) 

(9.1)  Determine  values  of  f^  (discussed  further  below.) 

(9.2)  Form  =  f^  by  (6.14). 


(9.3) 

(9.4) 

(9.5) 

(9.6) 


(9.7) 

(9.S) 

(9.9) 


(Tl  CT')  i  fT) 

Compute  d\  '  =  n.  ‘  a.  and  choose  d',^v 
J  1  j 

terminate  if  all  >  0. 


min  d 


(T) 

j 


<  0  or 


Compute  a 
Choose  9 


s(T)  "h 


^i(T-l)  ^h^^^  (6.11). 


r(T) 


by  (8.4)  or  terminate. 


Compute  by  (6.6,7),  transform  to  by 


one 


step  as  in  (6.11)  and  record  a(T)  for  j 


(T)  . 

r(T)  * 


(Optional)  Print  results  of  iteration  (See  Part  B) . 

i(  T) 

Condense  and  store  t)^''  and  r(T) . 

(Optional)  Check  solution  and  print  (See  Part  B) . 


X  -  THE  COMPOSITE  ALGORITHM  (Forming  f^, ) 

Suppose  that  a  basic  solution  has  been  obtained  which  is  infeasible, 
that  row  p  is  the  current  optimizing  form  and  that  q  is  the  smallest  row 
index  of  the  actual  restraint  equations. 


(10.1)  pj!; 


i(T-l)  „h(T-l) 


b"  ,  h 


CF  if  h  >  q  and  ^ 


Suppose  a  vector  a  /  v  has  somehow  been  chosen  to  bring  into  the  basis  and 

s(,  rj 

that  (8,4)  is  used  to  deterraj.ne  r(T).  After  the  change  of  basis,  the  new 
values  of  p^  are 

s(T). 


(10.2) 


j3^(T)  a,  0^^^^  >  0  (the  value  of  x°''^'') 


,i(T) 


(10.3)  (i  /  r(T)  ) 

Now  for  i  >  q  and  i^  F,  (10.3)  gives  >  0.  However,  for  i  £F,  there  are 


two  cases. 
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■*  i(T)  i(T-l) 

(10.4)  If  i€F  and  <  0,  then  p  ^  ^  >  P 

(10.5)  If  i£F  and  >  0,  then 

Hence  by  (10.2)  and  (10.3),  no  infeasibilities  are  created.  By  (10.2)  and 

(10.4),  some  infeasibilities  are  improved  or  removed  altogether.  However, 

by  (10.5)  some  may  be  made  worse.  To  overcome  this  difficulty,  consider  the 

function  X  =  ^  P^  <  0  which  is  a  measure  of  the  infeasibility  of  a  solu- 

i£  F 

tion.  We  wish  to  maximize  it  to  zero  and  hence  we  may  replace  the  maximiza¬ 
tion  of  p^  by  the  maximization  of 

(10.6)  o  -  p^'+  \  <  p^ 

provided  o  is  monotonically  non-decreasing.  When  o  reaches  its  maximm,  if 
X  =  0  then  p^  is  maximum.  If  4  <  0,  then  p^  is  too  great  and  X  must  be  in- 

.  p 

creased  without  regard  to  decreases  in  p  . 

Let  f.  =  1  if  i£  F  or  i  -  p,  f^  =  0  otherwise.  Then,  as  can  be  seen 

from  (9. 2, 3, 4) 


(10.7) 


4'^'  ■  ‘'i  4(t)  <  °  • 


Consequently,  after  change  of  basis,  the  new  value  of  o  is 

(10.8)  . 

If  the  set  F  is  void,  then  (10.8)  is  the  same  as 

(10.9)  -  V,) 

which  is  the  usual  formula  for  the  change  in  the  raaximand. 

Now,  however,  if  all  >  0,  we  cannot  claim  Terminal  2  without  check- 

ing  to  see  whether  F  is  void.  If  F  is  not  void  and  all  dj^''  >0,  it  may  be 
because  (10.9)  dominates  (10.8).®  In  this  case,  we  may  scale  down  fp  ,  perhaps 

*  Note  that  if  <  0,  then  <  0  for  i£F  and  i  ^  r(T). 


That  is,  ^  |fj^ 


to  zero,  until  such  time  as  \  =  0.  If  all  >  0  with  ”  0,  then  no 

feasible  solution  exists  since  X.  <  0  is  at  a  maximum. 

Even  though  X  =  0,  it  is  sometimes  desirable  to  set  f^  <  1,  Of  course, 

tolerances  must  be  built  into  the  code  in  testing  d.  >  0  since,  if  it  is 

sufficiently  small  in  magnitude,  it  ought  to  be  considered  zero.  Varying 

f  has  the  effect  of  varying  this  built-in  tolerance.  Provision  has  also 
P 

been  made  for  setting  f p  <  0  which  causes  to  be  minimized  instead  of  maxi¬ 
mized.  This  is  often  convenient  when  experimenting  with  a  model.  The  value 

of  f  is  entered  on  a  binarj''  card  subject  to  a  switch, 

P 

If  o  <  is  at  a  maximum  and  f  0,  the  machine  will  set  f  *  0  and 

P  P 

stop  so  that  other  values  may  be  loaded  if  desired.  If  X  “  0  and  f^  »  0, 

the  machine  will  set  f  =  1  and  stop. 

P 

The  variable  is  checked  for  monotonic  behavior  each  iteration,  accord¬ 
ing  as  f  >  0  or  f  <0.  This  test  is  suspended  if  X  <  0  or  when  making  arbi- 
P  P 

trary  transformations  when  the  behavior  of  p^'is  unpredictable. 


XI  -  PARAMETRIC  PROGRAMMING  (PLP) 

Provision  can  be  made  for  parametrizing  the  right  hand  side  as  a  linear 
function  of  0  >  0,  i.e. 

(11.1)  a^  x*^  *»  b^  +  0  c^ 

where,  if  a  Phase  I  was  used,  must  be  formed  in  the  same  manner  as  b^  in  (5 
To  do  this,  first  find  an  optimal  solution  for  0=0,  say 


Let 


(11.2)  ph(T-l)  ,  ^ 

(11.3)  c*"  =  o 


Now  using  Y^  place  of  ^  value  0^^^^  =  60  can  be  found  such 


.26). 


that 
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(11.4)  Y^  )  -  b^  +  C^ 


with 


(11,5)  >  0  for  i  >  q 

and 

p‘'(^-^>  -  Sr(T)  /  "  0  (r-r(T)). 

The  parameter  0  cannot  be  increased  by  more  than  with  the  basis 

^i(T-l)  without  violating  (11.5)  but  (11.4)  is  an  optimal  feasible  solution 

to  (11.1)  for  this  value  of  0.  Let 


(11.6) 


bUT)  .  hUT-1)  , 

1 

qKt)  ^  .KT-i)  _  e 


’r(T)  ^ 


TO  increase  0  further,  removed  from  the  basis  and  replaced 


i  1(T) 

with  some  a  to  form  a  new  basis 


(11.7) 


SO  that 


is  also  an  optimal  feasible  solution  to  (11.1)  for  the  same  0,  The  ^/hole 
process  can  then  be  repeated. 

The  index  s(T)  is  determined  by  the  ru.le  used  in  the  dual  simplex  algo¬ 
rithm.  Let  D  be  the  set  of  indices  j  for  which  ° 


Then  choose  <P  /'rp\  6y 


(11.8)  (j) 


P(T-I) 


s(T)  JiD  .  ,r(T-l)  h 


>0  (r  =  r(T)  ) 


Note  that 


and  (11,8)  is  the  analogue  in  the  dual  pi'oblem  of  (8.3). 

If  the  set  D  is  void,  then  0  is  at  a  maximum.  If  all  y  <0  (11.3), 

then  0  can  be  increased  without  bound.  These  are  the  only  two  automatic 
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terminations  in  PLPo  The  following  theorem  is  of  interest. 

Theorem;  If  the  choice  of  r(T)  for  (11.4)  is  t^nique  and  if  D  is  not  void, 

then  there  exists  a  finite  range  of  9,  ®  which 

the  solution  obtained  by  replacing  ^s(T)  *  ®(’^) 

determined  by  (11.8),  is  both  feasible  and  optimal, 

A  proof  can  be  found  in  Reference  4. 

A  separate  control  code  for  the  computer  is  used  in  doing  PLP.  It  always 
starts  from  a  prior  optimal,  feasible  solution.  Dub  to  (11.8),  the  iterations 
are  longer  than  the  regular  code. 

XII  -  MULTIPLE  PHASES 

It  is  somewhat  difficult  to  parametrize  the  optimizing  form  since  the 
analogue  of  would  be  a  row  vector  of  n+1  elements.  As  an  alternative, 
provision  is  made  for  multiple  optimizing  forms  which  can  be  made  to  differ 
by  finite  amounts  in  any  desired  way.  Of  course,  it  is  not  contemplated  that 
two  such  forms  will  be  drastically  different  since  that  is  equivalent  to  two 
different  problems.  In  such  a  case,  it  would  be  better  to  start  the  second 
one  from  the  beginning  or  from  the  end  of  Phase  I. 

It  is  also  possible  to  split  Phase  I  into  multiple  phases.  This  will 
not  be  discussed  further  since  its  application  is. limited  and  it  generalizes 
easily  from  the  discussions  given. 

It  will  be  easier  to  describe  the  use  of  multiple  phases  if  a  specific 
example  is  used.  Let  it  be  required  to  optimize  three  forms  and  to  start  the 
problem  with  a  Phase  I.  Then  the  auxiliary  form  must  be 

(12.1)  x'^  +  x^"^^  "0  (j“3,  4,  >...>n). 

The  form  to  be  optimized  first  (after  Phase  I)  must  be 

(12.2)  x^  +  a j  x*^  *>  0  (maximize  x^) . 


P-810 

3-14-56 

-27- 


Similarly,  the  other  two  forms  to  be  optimized  in  turn  must  be 

(12.3)  +  a j  «  0  (maximize  x^) 

(12.4)  x*^  +  a^  x*^  =«  0  (maximize  x*^)  , 

The  initial  restraint  equations  will  be 

(12.5)  x'^  +  »  b^  (i  -  4,  5»  mj  j  “  3,  4,  <>,,,  n). 

Thus  the  initial  value  of  q  is  specified  as  q  -  4  (total  number  of  phases) 
giving  p  *  q-1  =*  3  (number  of  "  Phase  II's"  ),  The  other  two  parameters 
required  are  z  ■■  1  (number  of  •'  Phase  I's")  and  ^  “  3  (index  of  sum  row,) 

At  the  end  of  Phase  I,  p,  q  and  z  will  all  be  reduced  by  1  giving 

2  10 
p«2,q  =  3>z**0  so  that  x  will  now  be  maocimlzed  disregarding  x  and  x  , 

(12,1)  will  now  be  considered  a  restraint  equation  the  same  as  (12,5),  and 

the  phase  will  be  recognized  as  a  Phase  II,* 

When  x'^  reaches  a  maximum,  p  will  be  reduced  by  1  to  p  =  1  but  q  will 

1  2  0 
remain  at  3  so  that  x  will  be  maximized  disregarding  x  and  x  ,  Similarly, 

when  x^  reaches  a  maximum,  p  will  be  reduced  to  p  ■=  0  with  q  still  remaining 

0  2  1 
at  3  so  that  x  will  be  maximized  disregarding  x  and  x  ,  In  other  words, 

p  is  reduced  each  phase  but  q  is  reduced  each  phase  only  as  long  as  z  can 

also  be  reduced.  All  three  must  remain  non-negative,  obviously. 


A  Phase  I  is  terminated  when  the  variable  being  maximized  reaches  zero, 

A  Phase  II  terminates  when  all  d^  >  0,  These  criteria  are  not  always  equi¬ 
valent  even  in  Phase  I,  If  is  representable  with  fewer  than  m  of  the 
columns  a^  (i,  j  >  0),  then  several  artificial  columns  may  remain  in  the 
basis  at  zero  level  at  the  end  of  Phase  I,  In  this  case,  the  d..  corresponding 

to  the  Phase  I  pricing  vector  will  not  all  be  non-negative.  If  Phase  II 

i  i 

starts  with  artificial  coliamns  in  the  basis  (other  than  a^  and  then 

a^_^^  may  be  eliminated  in  Phase  II  but  one  a^^^^  will  always  remain. 
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PART  B  -  0Ra4NIZATI0H  OF  THE  CODES 

We  will  describe  the  organization  of  the  present  programs  for  the  IBM 
704  which  are  a  culmination  of  the  experience  previously  discussed.  How- 
ever>  the  discussion  will  be  tied  to  the  machine  only  to  the  extent  of 
allowing  specific  statements  to  be  made.  Alternate  techniques  for  other 
machines  should  be  apparent  to  anyone  with  sufficient  familiairlty  with  comr- 
putors  to  have  a  real  interest.  The  equipment  assumed  available  is' the 
standard  704  with  4,096  words  of  core  storage,  four  drums  of  2,048  words 
each,  card  reader,  card  pvinch,  printer  and  five  magnetic  tape  units.  The 
data  assembly  program  will  be  discxissed  only  briefly o 
The  following  parameters  must  be  specified; 

m,  the  total  number  of  restraint  equations  including  all  but 
the  last  optimizing  form  (row  0.)  The  maximum  m  is  255. 
q,  the  total  number  of  phases, 
z,  the  niomber  of  Phase  I’s.  (Usually  0  or  1.) 

the  index  of  the  auxiliary  form  (sum  row)  for  Phase  I,  blank 
indicating  none . 

T,  the  ntimber  of  arbitrary  transformations. 

The  data  assembly  pi-ogram  reads  an  identification  card  followed  by  a  card 
with  these  parameters  and  then  by  cards  containing,  in  order: 

Initial  basis  headings,  if  any,  in  which, any  vectors  to  be  arbi¬ 
trarily  introduced  are  also  specified  by  making  negative. 

The  right  hand  side  b^. 

Auxiliary  right  hand  sides  c^,  if  any. 

The  elements  by  columns. 

Only  non-zero  elements  are  entered,  with  their  proper  indices.  The  program 
punches  in  binary  cards; 
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The  identification  card  (in  binary  coded  decimal-alphabetic) 

A  self  loading  origins  card  needed  by  the  loader  for  the  main 
routines  (2(m+l),  origins  for  H-,  V-,  V/-,  T-regions.  See 
below. ) 

That  part  of  K-region  which  is  pecxiliar  to  the  problem.  (See 
below. ) 

The  b^  in  double  precision,  floating  point . 

The  Sj  matrix  is  packed  and  written  on  tape  in  single  precision,  fixed  point, 
indexed  form.  It  may  also  be  pvinched  on  cards  if  desired. 

There  are,  of  course,  provisions  for  various  circumstances  and  certain 
built-in  checks,  as  for  example,  computing  and  inserting  the  sum  row  in 
a^  and  checking  tjnat  no  row  index  is  greater  than  m.  The  main  point  is  that 
considerable  thought  must  go  into  the  planning  of  this  program  to  mke  the 
operation  of  the  main  codes  as  efficient  as  possible.  We  turn  to  them  now. 

Conceptually,  the  high-speed  store  (HSS)  is  divided  into  two  main 
sections;  (i)  the  programs,  constants,  parauneters  and  temporary  storage, 

(ii)  the  data,  both  original  and  transformed. 

Similarly,  auxiliary  storage  is  divided  into  two  parts,  one  part  as  perman¬ 
ent  storage  for  (i),  the  other  as  permanent  storage  for  (ii).  Thus  the 
entire  HSS  is  in  a  sense  "temporary  storage"  .  In  practice,  a  considerable 
part  of  section  (i)  of  HSS  remains  intact  throughout  a  run  but  can  be  re¬ 
stored  from  auxiliary  storage  (drums)  if  necessary.  The  advent  of  the 
extremely  reliable  core-storage,  however,  makes  the  need  of  restoring  unlikely. 
When  HSS  reaches  the  sizes  now  contemplated  (or  as  much  as. .,12, 238  wortls), 
auxiliary  storage  for  the  code  and  for  certain  data  will  be  unnecessary.  The 

main  use  for  auxiliary  storage  is  for  the  matrix  and  the  vectors. 

j  r* 

Magnetic  tapes  are  used  for  both.  On  the  IBM  704,  three  tapes  are  actually 
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used  for  the  vectors  but  on  some  tape  systems  this  woxild  be  un¬ 


necessary. 

The  data  section  of  HSS  is  divided  into  four  regions: 
the  H-region  for  the  basis  headings;  m+1  words 
the  V-ragion  for  the  values  of  the  current  solution,  main¬ 
tained  in  double-precision,  floating  point;  2(m+l)  words 
the  W-region  for  work  space  in  generating  and  and  other 


purposes;  2(ni+l)  words. 

the  T-region  for  temporarily  holding  the  matrix  or  ^  vectors 
or  as  much  as  possible  of  either:  the  remainder  of  section 


(ii)  which  should  be  at  least  4(m+l)  words. 

The  size  of  these  regions  is  a  function  of  the  number  of  restraints.  Their 
origins  are  computed  by  the  data  assembly  program.  In  PLP,  V-region  is 
used  as  a  second  W-region  in  pricing.  A  duplicate  of  H-  and  V-reglons,  as 
well  as  the  original  b^,  is  kept  in  auxiliary  storage.  This  allows  re-starting 
after  an  error  and  checking  a  solution  by  computing  and  printing 


^1  .  ^l(T)  ph(T)  .  5,1. 

For  PLP,  it  is  also  necessary  to  keep  c^  and  in  auxiliary  storage. 

The  program  section  of  HSS  is  divided  into  seven  regions; 

(a)  Temporary  storage  called  COMMON. 

(b)  The  main  control  region,  called  the  MGR. 

(c)  A  sub- routine  for  doing  double  precision,  floating  point  addition, 
called  DPFADD. 

(d)  A  similar  sub— routine  for  multiplication,  called  DPFMUL. 

(e)  A  routine  called  the  distributor  (DISTRB) ,  explained  below. 

(f)  Space  for  the  largest  sub-routine.  The  first  location  is  called 


SRORIG, 
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(g)  Universal  constants,  parameters,  origins  etc.,  called  K-region. 

The  need  for  COMMON,  DPFADD,  DPFMUL  and  K-region  is  obvious.  Their  loca¬ 
tions  are  permanently  fixed  and  may  be  referred  to  from  any  program  in  the 
system.  All  sub-routines  are  closed.  The  K— region  contains  certain  cells 
whose  contents  are  fixed  only  for  part  of  an  iteration  or  other  sub-sequences 
of  the  piToblem.  There  are  conventions  on  when  and  by  what  routine  these  are 
to  be  changed.  COMMON  is  always  available  to  any  routine. 

The  function  of  the  MCR  is  to  make  the  major  decisions  and  call  for 
the  proper  sequence  of  operations.  Most  of  the  actual  work  is  delegated  "■ 
to  sub-routines.  Besides  the  MCR  for  the  main  algorithm,  there  is  one  for 
PLP  and  one  for  re-inverting  a  basis.  Other  auxiliary  programs  are  easily 
created  by  coding  a  new  MCR. 

The  MCR  calls  for  a  major  operation  by  linking  to  DISTRB  with  a  pseudo¬ 
operation.  DISTRB  calls  in  the  proper  sub-routine  from  auxiliary  storage 
and  loads  it  into  HSS  starting  at  SRORIG.  Control  is  then  transferred  to 
the  sub-routine  which  returns  control  to  the  MCR  after  completing  its  func¬ 
tion.  If  some  other  arrangement  for  linking  to  sub-routines  is  desired  (as 
for  example,  when  HSS  is  very  large)  it  is  only  necessary  to  change  DISTRB 
to  arrange  for  the  proper  routine  to  take  over  in  the  proper  way. 

There  are  fourteen  standard  sub-routines.  Others  may  be  added  for 
special  purposes  if  desired,  up  to  the  limit  of  auxiliary  storage  to  hold 
them.  More  importantly,  if  an  improved  or  a  special  version  of  one  is 
developed,  it  can  replace  the  old  one  merely  by  exchanging  the  proper  binary 
cai^s.  These  sub-routines,  together  with  DPFADD,  DPFMUL,  K-region  and  a 
special  loading  routine,  exist  in  binary  cards  which  constitute  a  basic  deck. 
The  origins  card  produced  by  the  data  assembly  is  put  ahead  of, and  the  MCR 
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behind  this  basic  deck,, 

As  soon  as  the  special  loading  routine  is  in  HSS,  it  takes  control 
and  loads  the  remainder  of  the  basic  deck  and  the  MCR>  storing  them  in 
auxiliary  storage  (drums).  Since  many  of  the  sub-routines  require  addresses 
to  be  set  which  depend  upon  m,  a  loading  interlude  is  perfonned  to  do  this 
initialization.  This  is  accomplished  by  an  initializing  routine  which  goes 
with  a  sub— routine^  the  two  being  loaded  simultaneously  into  HoS.  Control 
is  then  sent  to  the  initiaO-izing  routine  which  does  its  job  once  and  for  all 
and  returns  control  to  the  loader  .  The  loader  stores  the  initialized  sub¬ 
routine  in  auxiliary  storage  and  builds  a  catalogue  of  locations  into  DISTRB. 
The  MCR  takes  over  control  when  loading  is  complete. 

The  fourteen  sub— routines  (called  "  codes*)  are  as  follows.  Some  are 
used  for  multiple  purposes  which  are  controlled  by  internal  switches. 

Code  1.  Load  the  binary  cards  produced  by  the  data  assembly  and 

store  in  appropriate  places.  If  the  cards  are  included,  they 
are  transferred  to  tape.  Subject  to  a  switch,  a  value  may  be 
loaded  from  a  special  card  for  fp.  By  means  of  an  internal  switch 
and  the  contents  of  K— region  which  it  loads,  this  load  routine  is 
able  to  distinguish  various  cases,  initial  start,  re-start,  start 
of  PLP,  re-start  of  PLP,  so  that  the  proper  information  can  be 
stored  in  auxiliary  and  the  tapes  positioned  properly. 

Code  2.  Store  the  current  solution  in  auxiliary  (drums).  The 
•’  current  solution*'  is  defined  as  K-,  H-,  and  V-regions. 

Code  3.  (Normal)  Form  the  row  f^  in  W-region,  and  record  whether  any 
variables  are  infeasible. 

(During  PLP)  Form  sP  in  V-region  and  in  W-region. 

*In  so  far  as  the  linear  programming  codes  under  discussion  are  concerned, 
initisLlization  is  an  inovation  in  the  704  system  due  to  Hal  Judd  ox  lEM. 
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^  T  ^  i ( T— 1 ) 

Code  4.  (Normal)  Compute  V.'- region. 


(During  PLP)  Compute 


^(T-1)  v-region  and  ^^in  Vl-region. 


i  (T) 

Code  5.  (Normal)  Price  the  matrix  a.  and  choose  d^  . 

V 

(During  PLP)  Choose  ®-s  described  in  Part  A,  Sect.  XI. 

Subject  to  an  external  switch,  recognise  "  Curtain ••  marks  in  . 

Subject  to  instructions  from  the  MCR,  makd  the  following  check; 

(Normal)  If  d^"^)  /  0,  (during  PLP)  if  aj  /  0  for  j  / 

then  test  for  j  occurring  in  ^ 

occurred.  This  check  is  very  valuable  for  detecting  errors  before 
much  morfe  ccsnputing  is  done.  It  does  consume  some  amount  of  time, 
however,  which  increases  with  n. 

Code  6.  Load  a^^^^  from  aj  matrix  into  specified  region  as  double¬ 
precision,  floating  point  vector. 

Code  7.  Multiply  %(T)  vector  in  specified  region). 

Code  8.  (Normal)  Choose  the  index  r(T)  as  discussed  in  Part  A, 

Sect,  VIII.  For  arbitrary  transformations  or  inversion,  choose 


r(T)  by: 


|“s(T)|  "  ieA  |“s(T)| 
where  A  is  the  set  of  indices  i  for  which  (artificial). 

Code  9.  (Normal)  Compute  from  and  transform  to 

(PLP,  first  entriO  Compute  and  from  r(T)  and 

^i(T-l)^  (Part  A,  Elect.  XI).  Then  compute 


max  i 

•  -  _  or  . 
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i(T^  i 

(PLP,  second  entry)  Compute  '  from  transform 

^i(T-l)  ^l(T)^  (All  these  operations  are  very  similar,  even 

more  than  appears  at  first  glance.  The  variations  are  merely 
switches.) 

i(  T) 

Code  10.  Delete  zeros  from  n  '  ‘  and  index  non-zero  elements.  Store 

'r 

condensed  vector  in  auxiliary  storage  (dnims;  cf.  Code  13)* 

Code  11.  Multiply  out  taking  from 

i  f  T')  i 

original  a^  matrix  by  referring  to  ^  left  in  T-region 

and  b^  in  W-region  for  printing. 

Code  12,  Print  program.  This  program  is  quite  elaborate  and  longer 
than  the  other  codes.  Special  provisions  for  it  need  not  be  dis¬ 
cussed  here  except  to  say  that  appropriate  captions  and  identifi¬ 
cation  of  columns  are  printed  for  each  type  of  print-out,  of  which 
there  are  14,  The  print  output  can  also  be  put  on  the  fifth  tape, 
if  desired,  for  later  printing.  In  all  print-outs,  there  are  5 
columns  of  which  the  first  three  are;  and  i  «  0,  ..., 

The  last  two  are.  tt»  contents  of  W-  and  T-regions. 

Code  13.  This  performs  an  •%nd-of-stage  ”  in  which  some  of  the 
Kt") 

T)^'  '  are  transferred  from  drm  to  tape.  See  below. 

Code  14.-18.  (Undefined). 

Code  19.  Automatic  restart  program  for  recovering  after  an  error  by 

returning  to  the  beg.inning  of  the  iteration.  Clearly,  this  is 

highly  dependent  on  the  particular  machine.  The  important  points 

to  note  are  that  all  programs,  the  original  data,  the  current 

i(T) 

solution  at  beginnin.a:  of  the  iteration,  and  the  rip 


vectors  must 
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be  intact  somewhere  in  the  machine.  Also  this  code  must  have 
been  previously  instructed  where  to  return  control  to  and  be  able 
to  restore  all  eqviipment  (e.g.  tapes)  to  the  proper  positions. 

It  must,  of  course,  be  activated  by  some  external  means. 

Besides  the  physical  organization  above,  there  are  also  dynamic  groupings, 
several  of  which  have  already  been  discussed:  the  iteration,  the  phase, 
maximization  of  A<0,  arbitrary  transformations,  etc.  Another  grouping  is 
called  a  stage  and  always  consists  of  an  integral  number  of  iterations.  It 
has  nothing  to  do  with  phases  or  other  mathematical  aspects  but  is  simply  an 
operational  device.  It  was  devised  for  the- IBli  701  and  has  been  carried 
over  to  the  704,  for  slightly  different  reasons.  Some  analagous  arrangement 
is  probably  necessary  on  any  machine  with  non- homogeneous  storage  n^dia. 

Most  of  the  advantage  of  condensing  the  vectors  would  be  lost  if 

a  separate  access  to  auxiliary  storage  had  to  be  made  for  each  one.  It  is 
desirable  that  as  big  a  chunk  of  these  transformation  vectors  as  possible  be 
recalled  at  one  time.  The  limiting  factor  is  the  size  of  T-region.  Hence, 
as  the  are  generated,  they  are  stored  on  a  drum  until  one  more  would 

exceed  the  capacity  of  T-region.  At  this  point  an  end-of-stage  procedure  is 
performed.  Before  describing  this,  it  is  necessary  to  explain  the  use  of 
three  tapes  for  the  vectors. 

The  tapes  on  the  704  can  be  back-spaced  and  additional  records  can  be 
added  to  an  existing  file.  However,  they  cannot  be  read  in  a  backward 
direction  and  the  back-space  is  somewhat  slow.  Hence  the  records  are  stored 
in  reverse  order  on  a  second  tape  for  use  in  ccnnputing  •  The  third  tape 
is  used  alternately  with  the  second  from  stage  to  stage.  Though  this  costs 
some  copying  time,  it  does  provide  a  spare  (and  checked)  copy  of  the  tape  at 
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last  st»ge  but  one  which  is  very  desirable.  It  is  dangerous  for  one  to  burn 
all  his  bridges  behind  him,  especially  with  magnetic  tapes. 

The  end-of-stage  procedure  is  as  follows: 


(1) 

(2) 

(3) 


Compute 

Print  i,  b^,f 


Transfer  the  on  drum  to  an  additional  record  on  the 

forward  tape. 


(4)  Transfer  the  on  drum  to  the  first  record  on  tte  free 

backward  tape. 

(5)  copy  old  backward  tape  to  remainder  of  new  backward  tape. 

(6)  Punch  out  binary  cards  containing  K-region,  and 

During  PLP,  also  p\mch  and 

.  ,  (7)  Adjust  the  necessary  bookkeeping  parameters o 
An  end-of-stage  may  be  forced  by  an  external  switch.  It  also  occurs  at  the 
end  of  Phase  I  6ind  terminations.  Steps  (1)  and  (2)  above  may  be  forced 
without  the  others  at  the  end  of  a  cycle  by  an  external  switch.  This  would 
be  desirable  on  machines  which  do  not  require  the  above  tape  manipulations. 

To  restart  from  the  end  of  a  stage,  it  is  only  necessary  to  use  the 
punch-outs  to  replace  the  corresponding  original  data  cards  and  put  the 
tapes  back  on  the  same  units. 

Similar  simple  hand  collating  of  blocks  of  punch-outs  with  the  original 
data  cards  (binary)  and  changing  of  the  MCE  are  all  that  are  necessary  to 
set  up  the  deck  for  PLP  or  re-inverting  the  basis.  .Special,  short  MCR's  are 
usually  accumulated  for  such  things  as:  transforming  a  new  right  hand  side, 
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computing  and  printing  the  d^,  altering  a  vector  in  the  basis,  etc.  A 

routin,  13  provided  for  copying  a  tap.  In  reverae  order  ao 

that  a  problea  can  be  picked  up  after  bad  luck  with  tapes. 
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PART  C  -  PRS-DESIGWED  MODELS  AND  ALTERNATE  METHODS 


The  objectives  which  have  motivated  the  programs  discussed  in  Parts  A 
and  B  have  been  mainly  three  fold; 

(i)  to  perform  calculations  for  any  linear  programming  model 

which  might  be  presented,  wi.thin  current  computing  capacity, 
in  as  efficient  and  accurate  a  way  as  possible; 

(ii)  to  accept  any  reasonable  special  requirements  and  even  to 
anticipate  them; 

(iii)  to  increase  the  size  of  problems  which  are  computationally 
feasible . 

Even  assuming  a  degree  of  success  in  all  three,  there  is  much  left  to  be 
desired.  As  to  (ii),  it  is  obvious  thfit  new  requests  will  continually  be 
made  which  are  unexpected.  V/hile  a  few  special  devices  can  be  pre-designed, 
the  main  reliance  must  be  on  an  extremely  flexible  code  in  which  the  algorithm 
can  be  modified  in  any  particular  without  upsetting  the  whole  apple-cart. 

We  believe  we  are  in  a  better  position  in  this  regard  than  in  the  other  matters. 
Although  it  takes  time  to  carefully  plan  the  organization  of  a  computer  pro¬ 
gram,  it  pays  off  in  the  long  run.  The  important  thing  is  that  the  code 

truly  reflect  the  basic  method  and  not.be  a  hodge-podge  of  standard  library 

sub-routines  assembled  for  expediency.  This  will  be  even  more  important  for 
the  higher  levels  of  abstraction  which  will  be  required  in  some  proposed 

methods.  It  is  also  important  for  the  programmer  t-o  be  on  the  look-out  for 

situations  amenable  to  simple  tricks.  For  example,  dualizing  a  problem 
sometimes  renders  special  computations  easier.  Incidentally,  we  have  not 
found  the  actual  dual  algorithm  very  useful  in  practice.  It  is  slow  and  the 
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compoaite  algorithm  was  developed  mainly  to  replace  its  use.  To  a  lesser 
extent,  the  same  is  true  of  PLP. 

It  seems  that  (i)  and  (iii)  are  somewhat  antithetical.  Accuracy  is 
usually  not  an  important  consideration  in  final  restilts  but  is  necessary  to 
maintain  the  progress  of  the  computations.  Efficiency  is  a  difficult  thing 
to  measure  and  depends  on  one's  perspective,  but  for  very  large  problems  it 
can  be  taken  as  meaning  lowest  possible  machine  time.  Hence,  accuracy  and 
efficiency  are  as  much  a  part  of  (iii)  as  of  (i).  But,  if  the  size  of 
problems  is  to  be  increased  efficiently,  then  special  methods  must  be  developed 
and  most  of  the  suggestions  so  far  depend  on  structure  in  the  a^  matrix. 

t) 

Thus  it  appears  that  we  must  begin  to  insist  that,  if  problems  are  to  exceed 
a  certain  sice  (and  this  has  been  the  main  demand),  then  they  must  fit  into 
certain  patterns.  This  my  not  be  as  restrictive  as  it  sounds.  It  seems  un¬ 
likely  that  many  problems  would  be  formulated  by  someone  just  writing  dovm, 
say,  500  random  equations  as  a  set  of  constraints.  Indeed,  we  anticipate 
that,  as  systems  increase  in  size,  there  will  be  more  requests  for  constructing 

the  matrix  itself  by  machine.  If  tiiis  is  so,  then  it  virtually  implies 
d 

a  pattern  and  there  should  be  some  degree  of  freedom  in  the  way  it  is  to  be 
constructed. 


The  nv3Lin  difficulty  at  present,  of  course,  is  to  devise  a  set  of  patterns 
which  lead  to  efficient  computing  procedures  and  which,  at  the  same,  time  en¬ 
compass  realistic  problems.  Work  along  these  lines  has  barely  begun.  :le  will 


♦Markowitz  has  a  straight-forv/ard  method  of  obtaining  a  new  feasible  solution 
after  a  change  in  the  right  hand  side,  provided  the  change  is  proportional  in 
a  positive  sense  to  legitimate  vectors  in  the  system.  This  has  not  been  used 
extensively  but  has  much  merit.  It  is  simply  the  application  of  the  theorem  on 
basic  solutions;  If  any  feasible  solution  exists,  then  a  basic  feasible  solu¬ 
tion  exists. 
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mention,  however,  one  such  pre-dasigned  model  which  is  now  under  development. 

In  numerous  applications  of  linear  programming,  it  has  been  observed 
that  the  matrix  can  be  partitioned  into  e.  special  block  triangular  form. 

•J 

Extending  our  notation  to  elements  which  are  thanselves  matrices,  let 


(1  1#  .a*,  >  J  *  ^9  •••# 

Then  the  structured  matrix 


A°  A°  A^ 
Hq  ^2 


-  ^ 


(unspecified  blocks  all  zero) 


gives  rise  to  bases  which,  by  a  relatively  simple  transformation,  can  be 
maintained  in  a  similar  structured  form.  It  is  fairly  clear  that  most  of 
the  numerical  operations  can  be  confined  to  the  individual  blocks,  assuming 
the  corresponding  diagonal  blocks  of  the  basis  have  rank  m^,  respectively. 

A  groat  deal  of  study  has  gone  into  devisihg  the  proper  sequence  of  opera¬ 
tions  to  make  the  arithmetic  as  short  as  possible. 

Examples  of  this  structure  occur  in  the  transportation  problem,  the 
metal-working  problem  studied  by  Markowitz,  the  gasoline  blending  model  of 
Chames ,  Cooper  and  Mellon,  and  the  two-period  air-transport  procurement 
problem  of  Manne.  A  special  method  based  on  this  structure  would  not  be 
economical  in  all  instances  but  it  is  believed  there  is  considerable  potential 
demand  for  such  a  technique. 


P-810 

3-14-56 

-41- 

We  cannot  close  without  outlining  another  technique  which  has  recently 
been  used  with  spectacular  results  on  small  problems.  Like  the  one  above, 
it  was  proposed  in  an  attempt  to  increase  the  size  of  problems  which  can  be 
handled.  In  both  cases,  it  turns  out  that  these  methods  are  more  efficient 
for  some  pio>blems  well  within  the  capacity  of  current  codes  than  the  stand¬ 
ard  methods.  In  the  following,  there  is  no  special  assumption  on  the  a^ 

matrix  except  that  it  be  sparse.  The  necessary  degree  of  sparseness  for 
efficiency  is  not  yet  known  since  it  is  heavily  dependent  on  code  logic. 

In  1954,  Markowitz  pointed  out  that  if  it  were  possible  to  obtain  a 
form  of  inverse  v/hich  could  be  applied  qtiickly  and  which  was  not  too  expensive 
to  generate,  one  could  invert  a  basis  and  then  use  the  product  form  to 
modify  it  for  several  iterations.  When  the  transformations  became  too  long, 
then  the  latest  basis  could  be  re-inverted,  there  being  some  optimal  point 
for  doing  this. 

Now  it  is  the  number  of  multiplications  and  additions  that  determine 
the  time  for  applying  an  inverse,  in  short  the  number  of  non-zero  elements  in 
the  form  of  inverse  used.  It  is  theoretically  possible  to  obtain  a  form  of 
inverse  in  which  the  number  of  non-zero  elements  is  very  little  greater  than 
the  number  of  non-zero  elements  in  the  original  matrix.  The  product  form 
has,  in  general,  fewer  non— zero  elements  than  the  full  inverse  of  a  sparse 
matrix.  For  example; 


10  0  1 

-1 

1/2  -1/2  -1/2  -1/2 

-110  0 

a 

1/2  1/2  -1/2  -1/2 
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0 
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1^ 
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with  16  non-zero  elements.  However  A"^  can  be  represented  in  product  form 
by  the  fourV]^^^^  vectors,  where  r(t)  -  t: 


which  have  10  non-zero  elements.  Note  that  if  r(l)  “  4,  the  number  of 
elements  in  the  product  fom  changes.  Thus,  the  product  fom  does  not,  in 
general,  minimize  the  number  of  non-zero  elements  required,  in  fact  there  is 
no  known  algorithm  for  doing  so,  at  least  in  a  practical  sense.  However,  as 
Markowitz  further  pointed  out,  one  does,  in  fact,  do  a  pretty  good  job  of  it 
in  solving  a  system  of  simultaneous  equations  by  hand.  If  several  co¬ 
efficients  in  the  equations  are  zero,  one  does  not  proceed  to  triangilLarize 
in  a  straightforward  manner  but  rather  seeks  to  eliminate  first  those  rows  and 
columns  with  the  fewest  elements . 

His  first  proposal  was  to  operate  first  on  a  grid  or  "picture"  of  the 
matrix,  following  through  the  effect  of  the  eliminations  on  the  zeros,  with¬ 
out  actually  obtaining  the  numerical  values  of  the  transfomed  elements 

The  result  of  this  was  to  be  an  "  agenda"  consisting  of  a  sequence  of 

h(t)’ 

pairs  of  indices  (r(t),  s(t)),  t  ^  1,  2,  ...,  The  rule  for  choosxng 

r(t)  and  s (t)  was  to  be  as  follows: 

Let  be  the  count  of  the  non-zero  elements  remaining  in  column 

h  at  this  point  of  the  operations  and  let  be  similarly  the 

count  for  row  i.  Then,  r  and  s  are  chosen  by 

min  ‘ 

“h(t)  ^  °  ^ 
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Having  detenained  the  agenda,  the  actual  transformations  were  to  be  carried 

through  forming  ^forward  solutions "  and  "  backward  solutions "  in  the  usual 

i(t ) 

manner  of  eliminating.  The  forward  transformations  are  like  our 
matrices  and  the  backward  transfonnations  like  the  transpose  of  these. 

There  were  several  technical  difficulties  with  coding  up  such  a  pro¬ 
cedure  that  would  pi-oduce  a  result  compatible  with  the  product  form,  and 
the  bookkeeping  was  unusually  difficxilt.  There  was  also  a  theoretical 
difficulty.  Suppose  that  in  carrying  out  the  agenda,  one  of  the  pivot 

elements  became  zero  by  chance.  V.'hat  course  of  action  should  be  taken 

s(  t; 

and  could  anything  be  salvaged? 

The  project  was  shelved  for  several  months  and  then  re-investigated. 

It  was  found  that  the  agenda  and  the  actual  transfonnations  could  be  per¬ 
formed  simultaneously,  O'  rather,  in  alternating  steps.  If  a  pivot  element 
turns  out  zero  it  is  simply  rejected  and  another  one  picked.  Also,  by 
slightly  changing  the  procedure,  the  backward  as  well  as  the  forward  trans¬ 
formations  can  be  put  into  identically  the  same  form  as  the  usual  product 
form.  There  is  a  permutation  occurring  between  the  two  but  this  is  recorded 
simply  by  insartJjig  the  agenda  at  this  point.  The  complete  details  must 
await  a  later  paper. 

A  program  for  li/iND's  own  macldne,  the  JOHI^MIAC,  has  now  been  in  opera¬ 
tion  for  a  short  time.  It  has  given  remarkably  good  results  on  problems  up 
to  51  eqxmtions  in  some  kOO  variables.  The  code  is  limited  to  12?  equations 
but  this  figure  is  probably  ambikous  since  the  machine  has  no  tapes  at  present. 
There  are  4096  words  of  HSS  and  12,2m  words  of- drum  storage.  The  inversion 
code  works  so  well  that  we  have  designed  the  JOHUNIAC  simplex  code  so  that, 
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if  a  job  is  interrupted,  the  basis  must  be  re-inverted  before  continuing. 
Although  the  machine  is  slower  than  the  701,  this  code  is  faster  than  the 
701  codes,  partly  due  to  fixed  point  arithmetic  (double-precision)  and  very 
short  drum  transfer t imes.  However,  on  one  problem,  this  form  of  inverse 
had  only  107  per  cent  as  many  non-trivial  elements  as  the  given  basis.  This 
makes  the  inverse  fast  to  apply,  including  the  part  already  generated  vdien 
generating  the  inverse  itself,  that  is,  even  the  generation  of  the  inverse 
is  speeded  up  by  advantageous  feedback  of  its  own  features. 

It  is  too  early  to  say  how  successful  this  scheme  would  be  on,  say,  the 
704  with  a  very  large,  but  sparse,  system,  say  400  or  500  equations.  It  has 
been  amazing  to  the  writer  that  the  tremendous  amount  of  bookkeeping  and 
"  juggling''  in  the  JOHNNIAC  code  actually  pays  off.  However,  this  is  in  line 
with  all  our  other  experience,  that  one  can  afford  to  do  aliQOst  anyttiing  to 
cut  down  the  number  of  non-zero  multipliers  and  addends  recorded.  There  are 
still  one  or  two  minor  points  of  difficulty  to  be  cleared  up  in  the  method  but 
the  potential  flexibility  exceeds  even  that  of  the  704  codes  in  some  respects. 
(The  present  JOHNNIAC  codes  are  not  in  a  polished  form  but  are  quite  con¬ 
venient  to  use,  even  so.) 

There  are  thus  at  least  four  variations  of  the  simplex  method  for  gen¬ 
eral  problems,  any  one  of  which  might  be  most  advantageous  depending  on  the 
machine  and  the  size  of  problems  to  be  provided  for. 

The  original  simplex  method;  for  fairly  small  machines  (2000 
words  total  capacity)  and  small  problems,  this  uses  a 
fixed  amount  of  storage  and  a  relatively  simple  code. 

The  revised  simplex  method  with  explicit  inverse:  for  small 


machines,  the  size  of  problem  ch.n  be  increased  somewiiat 
by  keeping  only  the  inverse  in  HSS  and  muining  the 
matrix  through  on  cards,  for  example,  again  and  again. 

The  code  can  be  kept  fairly  simple  and  much  flexibility 
attained. 

The  revised  simplex  method  with  product  form  of  inverse;  for 
elaborate,  automatic  set-ups,  this  method  is  serving  v/ell 
but  has  definite  limitations  for  problem  size,  especially 
without  tapes. 

The  revised  method  with  product  form  combined  with  elimination 
■  form  of  inverse;  extends  the  size  of  piNiblems  that  can 
be  run  on  a  given  raaciiine  although  not  suitable  for  small 
machines  due  to  largo  amount  of  code.  The  main  drawback  is 
'  that  its  advantages  do  not  take  effect  until  a  ,non~trivial 

basis  is  knovnio 

It  is  probably  a  waste  of  time  to  attempt  to  run  problems  on  a  machine 
with  less  than  2000  words  of  storage.  No  matter  what  method  is  used,  the 
necessary  code  will  cut  down  available  storage  drastically  and  any  attempt  at 
packing  information,  or  other  such  tricks,  will  onlj'  make  the  code  longer  at 
the  further  expense  of  data  storage.  A  program  for  15  or  20  order  systems 
might  be  used  for  a  while,  but  before  long,  the  customers  would  v/ant  to  run 
bigger  problems. 

We  believe  that  the  time  has  come  for  specialization  in  the  field  and 
that  future  efforts  should  be  directed  toward  concentrated  study  of  problems 
of  particular  types,  both  coraputationally  and  an.  ..^v-tically.  If  the  experienc 
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and  attempts  related  in  the  foregoing  help  in  any  way  to  avoid  future  pit- 
falls,  they  will  have  served  one  of  their  purposes. 
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